R Notebook

library(tidyverse)
library(dplyr)
library(readxl)
library(broom)
library(MuMIn)
library(lme4)

Pred <- read_csv("../data/predator.csv.gz",na = c("","n/a","NA"),col_types = cols(`SD PP`=col_double()))
problems(Pred)

colnames(Pred) <- gsub(" ", "_", colnames(Pred))

Pred 

Pred$Latitude <- Pred$Latitude %>% 
    sub("º", "d",.)
library(sp)
Pred$Latitude <- as.numeric(char2dms(Pred$Latitude))

Pred <- Pred %>%
mutate(Specific_habitat=gsub("Coastal bay", "Coastal Bay", Specific_habitat)) %>%
mutate(Specific_habitat=gsub("shelf", "Shelf", Specific_habitat))

pred_abundance dataset

This dataset is a subset of the marine dataset. What does it include?:

# pred_abundance <- Pred %>% 
#     group_by(Predator) %>% 
#     tally() %>% 
#     arrange(desc(n)) %>% 
#     filter(n > 800)
# pred_abundance
# 
# prey_abundance <- Pred %>% 
#    group_by(Prey) %>% 
#     tally() %>% 
#     arrange(desc(n)) %>% 
#     filter(n > 50 & n < 3000)
# prey_abundance
pred_prey_dataset <- Pred %>%
    # filter(Predator %in% pred_abundance$Predator) %>%
    # filter(Prey %in% prey_abundance$Prey) %>%
    # filter(Depth < 1000) %>%
    group_by(Predator_common_name, 
             Prey_common_name, 
             Specific_habitat, 
             Type_of_feeding_interaction, 
             Predator_lifestage, 
             Mean_PP,
             Mean_annual_temp,
             Depth) #%>%
   # summarise(Latitude = mean(Latitude),
   #           Predator_mass = mean(Predator_mass),
   #           Prey_mass = mean(Prey_mass),
   #           Depth = mean(Depth))
     
  
pred_prey_dataset

log predator mass vs log prey mass

As prey mass increases so does predator mass. When plotted this way, the data clusters into 4-5 blobs – it’s possible that there might be something underlying this pattern.

pred_prey_dataset %>% 
    ggplot(aes(x = log(Prey_mass), y = log(Predator_mass))) +
    geom_point() +
    geom_smooth(method = "glm") 

NA

Colour by specific habitat

Same plot, but now we colour data based on what habitat the predator-prey interaction is. It appears that habitat predicts the observed pattern… so what is it about habitat that underlies this pattern that as prey size increases so does predator size?

pred_prey_dataset %>% 
    ggplot(aes(x = log(Prey_mass), y = log(Predator_mass), colour = Specific_habitat)) +
    geom_point() 

# based on this graph, prey mass also appears to be smaller than predator mass for most of the time (1:1 line, maybe get a stat on this)

Colour by feeding interaction

Same plot, only we colour the data by the kind of feeding pattern the predator exhibits. It looks like this could also explain the data, although maybe not as strongly as habitat – you see that predacious feeding spans the extremes, and picsivorous feeding does too. But these are really general feeding types… are feeding types found in more habitats than in others?

pred_prey_dataset %>% 
    ggplot(aes(x = log(Predator_mass), y = log(Prey_mass), colour = Type_of_feeding_interaction)) +
    geom_point() 

Predacious data

predacious <- pred_prey_dataset %>% 
    filter(Type_of_feeding_interaction == 'predacious') %>% 
    select(Predator_mass, Prey_mass, Type_of_feeding_interaction, Prey_common_name, Specific_habitat) 

Adding missing grouping variables: `Predator_common_name`, `Predator_lifestage`, `Mean_PP`, `Mean_annual_temp`, `Depth`

   
predacious %>% 
    ggplot(aes(x = log(Prey_mass), y = log(Predator_mass), colour = Prey_common_name)) +
    geom_point()

# looking at the predacious data, there's two extremes:
# squid prey seem to make predator prey really big
# on the other hand, copepods and invertebrates appear to make predators really small
# an energy transfer issue? 
predacious %>% 
    ggplot(aes(x = log(Prey_mass), y = log(Predator_mass), colour = Specific_habitat)) +
    geom_point() 

# this might just be a function of the ice zone having a lot of observations, but being in the ice zone or nearshore makes you small, while living in the shelf makes you large
# in general, the slope of this looks almost the same as the previous graphs 

Planktivorous data

planktivorous <- pred_prey_dataset %>% 
    filter(Type_of_feeding_interaction == 'planktivorous') %>% 
    select(Predator_mass, Prey_mass, Type_of_feeding_interaction, Prey_common_name, Specific_habitat)

Adding missing grouping variables: `Predator_common_name`, `Predator_lifestage`, `Mean_PP`, `Mean_annual_temp`, `Depth`

planktivorous %>% 
    ggplot(aes(x = log(Prey_mass), y = log(Predator_mass), colour = Prey_common_name)) +
    geom_point() 

planktivorous %>% 
    ggplot(aes(x = log(Prey_mass), y = log(Predator_mass), colour = Specific_habitat)) +
    geom_point()

# all planktivorous species are found in the ice zone, and they're all really small
# slope increases quickly 

Piscivorous data

piscivorous <- pred_prey_dataset %>% 
    filter(Type_of_feeding_interaction == 'piscivorous') %>% 
    select(Predator_mass, Prey_mass, Type_of_feeding_interaction, Prey_common_name, Specific_habitat) 

Adding missing grouping variables: `Predator_common_name`, `Predator_lifestage`, `Mean_PP`, `Mean_annual_temp`, `Depth`

piscivorous %>% 
    ggplot(aes(x = log(Prey_mass), y = log(Predator_mass), colour = Prey_common_name)) +
    geom_point()

# You see some grouping here, it's a little more spaced out than what you saw for predacious, but there is definitely some grouping
piscivorous%>% 
    ggplot(aes(x = log(Prey_mass), y = log(Predator_mass), colour = Specific_habitat)) +
    geom_point()

# Specific habitat is also a really good indicator for prey size - predator size

Facet wrapping

pred_prey_dataset %>% 
    ggplot(aes(x = log(Prey_mass), y = log(Predator_mass), colour = "gray50")) +
    facet_wrap(~ Type_of_feeding_interaction+Specific_habitat ) +
    geom_point() +
    geom_smooth(method = 'glm', colour = 'black')

#ggsave("facetwrapping.jpg", width = 50, height = 40, units = "cm")

Some latitude stuff

Pred %>% 
    ggplot(aes(x = log(Prey_mass), y = log(Predator_mass), colour = abs(Latitude))) +
    geom_point(alpha = 0.3)

pred_prey_dataset %>% 
    ggplot(aes(x = log(Prey_mass), y = log(Predator_mass), colour = Predator_common_name)) +
    facet_grid(Specific_habitat ~ Type_of_feeding_interaction) +
    geom_point() +
    geom_smooth(method = 'glm', colour = 'black')

# ggsave("facetwrapping_adults.jpg", width = 50, height = 40, units = "cm")

pred_prey_dataset %>% 
    filter(Type_of_feeding_interaction == "planktivorous") %>% 
    ggplot(aes(x = log(Prey_mass), y = log(Predator_mass))) +
    facet_wrap(~ Specific_habitat) +
    geom_point() +
    geom_smooth(method = 'glm', colour = 'black')

# ggsave("facetwrappingplanktivorous.jpg", width = 50, height = 40, units = "cm")

pred_prey_dataset %>% 
    filter(Type_of_feeding_interaction == "piscivorous") %>% 
    ggplot(aes(x = log(Prey_mass), y = log(Predator_mass))) +
    facet_wrap(~ Specific_habitat) +
    geom_point() +
    geom_smooth(method = 'glm', colour = 'black')

# ggsave("facetwrappingpiscivorous.jpg", width = 50, height = 40, units = "cm")

pred_prey_dataset %>% 
    filter(Type_of_feeding_interaction == "predacious") %>% 
    ggplot(aes(x = log(Prey_mass), y = log(Predator_mass))) +
    facet_wrap(~ Specific_habitat) +
    geom_point() +
    geom_smooth(method = 'glm', colour = 'black')

# ggsave("facetwrappingpredacious.jpg", width = 50, height = 40, units = "cm")

Notes from meeting

Before we were thinking Prey and Predator have a really strong correlation but there’s nothing really here! (I think it’s moer like, depending on where you are, the correlation might be strong but not necessarily strong. But if we have to take some habitats out, it’d be interesting to hypothesize why the slopes of some these things change at different rates – could it be temperature? could it be depth?). Predator is strongly correlated with their habitat, which has its own depth. WHAT YOU EAT, AND WHERE YOU EAT affects the correlation between predator and prey mass. GLM: look at interactions between specific habitat, type of feeding interactions?? Only when we looked at generalist eaters did habitat

Predator-prey body mass ratios

# pred_prey_dataset %>%
#     mutate(PPMR = mean(Predator_mass/Prey_mass)) %>% 
#     group_by(Specific_habitat, PPMR) %>% 
#     select() %>% 
#     ggplot(aes(x = Specific_habitat, y = log(PPMR))) +
#     geom_col() +
#     theme(axis.text.x = element_text(angle=-45, hjust=0, vjust=1)) 
# Not really sure if this tells us anything -- I think it gives us the same information as the slope 
# ggsave("predator prey ratio.jpg", width = 15, height = 10, units = "cm")

OKAY PROJECT

Question: How do specific habitats affect predator-prey body size relationships? (specifically, how does it affect log prey mass (dependent variable)) - In other words, when log predator mass and log prey mass are plotted, what explains the variation in strength of the slope? - Is there an interaction effect with feeding type? With predator mass?

Predator-prey relationships

Same graphs again:

pred_prey_dataset %>% 
    ggplot(aes(x = log(Predator_mass), y = log(Prey_mass), colour = Specific_habitat)) +
    geom_point(size = 0.3) 

pred_prey_dataset %>% 
    ggplot(aes(x = log(Predator_mass), y = log(Prey_mass), colour = Type_of_feeding_interaction)) +
    geom_point(size = 0.3) 

# Prey mass = dependent variable
# Predator mass = independent variable: Brose et al. 2010 use Predator mass as the independent variable because there is less error in it than prey mass (prey are usually removed from predator guts, sometimes hard to identify exact species)
# are there other independent variables? do they interact?
# Question: how do you add a 1:1 line?
# fit some models to this graph: mixed effects, AIC

pred_prey_dataset %>% 
    group_by(Type_of_feeding_interaction) %>% 
    tally()

pred_prey_dataset %>% 
    group_by(Specific_habitat) %>% 
    tally()

# I don't think we need to worry about the lifestages of predators becauase in the dataset it distinguishes between what the predator eats depending on its lifestage

filtered_dataset <- pred_prey_dataset %>% 
    filter(Type_of_feeding_interaction != 'insectivorous' & Type_of_feeding_interaction != 'predacious/piscivorous') %>% 
    filter(Specific_habitat != 'Coastal, SW & SE Greenland' & Specific_habitat != 'inshore' & Specific_habitat != 'Nearshore waters')
filtered_dataset 

pred_prey_dataset %>% 
    ggplot(aes(x = log(Predator_mass), y = log(Prey_mass), colour = Geographic_location)) +
    geom_point(size = 0.1) +
    facet_wrap(~ Specific_habitat) +
    geom_smooth(method = 'glm', colour = 'black')

# I think the slope might tell you about the range of body sizes, as well as how prey body mass changes with predator body mass in different habitats

pred_prey_dataset %>%
    ggplot(aes(x = log(Predator_mass), y = log(Prey_mass), colour = Type_of_feeding_interaction)) +
    geom_point(size = 0.1) +
    facet_wrap(~ Type_of_feeding_interaction) +
    geom_smooth(method = 'glm', colour = 'black')

Model

library(lme4)

filtered_dataset <- Pred

Linear model 1

This model says: That for every habitat, predator and prey mass relationships have a different intercept (different starting point), and the effect is has on the slope of the relationship is also different in every habitat. Different slope and intercept for habitat.

lm1 <- lmer(log(Prey_mass) ~ 
                log(Predator_mass) + 
                (1 + log(Predator_mass)|Specific_habitat), 
            data = filtered_dataset)
summary(lm1)

Linear mixed model fit by REML ['lmerMod']
Formula: log(Prey_mass) ~ log(Predator_mass) + (1 + log(Predator_mass) |      Specific_habitat)
   Data: filtered_dataset

REML criterion at convergence: 153073.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-6.3008 -0.5303 -0.0086  0.6094  7.2546 

Random effects:
 Groups           Name               Variance Std.Dev. Corr 
 Specific_habitat (Intercept)        19.3577  4.3997        
                  log(Predator_mass)  0.3522  0.5934   -0.09
 Residual                             4.6589  2.1585        
Number of obs: 34931, groups:  Specific_habitat, 16

Fixed effects:
                   Estimate Std. Error t value
(Intercept)         -5.6893     1.1100  -5.125
log(Predator_mass)   0.9350     0.1564   5.977

Correlation of Fixed Effects:
            (Intr)
lg(Prdtr_m) -0.073

ranef(lm1)

$Specific_habitat
                                (Intercept) log(Predator_mass)
Coastal Bay                        2.341380        -0.07530835
Coastal, SW & SE Greenland         3.972445        -0.52946851
demersal food web                  3.055145        -0.12506006
estuary/coastal                   -1.296394         0.09789005
Euboikos and Pagassitikos Gulfs    3.326412         0.16178262
inshore                           -3.834071        -0.24624064
nearshore                         -3.566458         0.07125153
Nearshore waters                  -1.766406         0.33218159
Off the continental Shelf         -8.556339         0.61526105
offShelf and on Shelf              5.862114         1.56673981
open ocean                         6.563611        -0.97400337
Pelagic                            2.143563        -0.22433598
Seasonal Pack Ice Zone            -6.103949        -0.40220529
Shelf                              2.662171        -0.44925209
Shelfbreak/open ocean             -2.623919         0.15605431
transition region                 -2.179304         0.02471333

tidy(lm1,conf.int = T)

Linear model 2

Assuming different slopes and intercepts for different habitats and feeding types.

lm2 <- lmer(log(Prey_mass) ~ log(Predator_mass) + (1 + log(Predator_mass)|Specific_habitat) + (1 + log(Predator_mass)|Type_of_feeding_interaction), data = filtered_dataset)

Model failed to converge with max|grad| = 0.00435131 (tol = 0.002, component 1)

summary(lm2)

Linear mixed model fit by REML ['lmerMod']
Formula: log(Prey_mass) ~ log(Predator_mass) + (1 + log(Predator_mass) |  
    Specific_habitat) + (1 + log(Predator_mass) | Type_of_feeding_interaction)
   Data: filtered_dataset

REML criterion at convergence: 147365.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-6.3194 -0.5476 -0.0219  0.5822  7.7991 

Random effects:
 Groups                      Name               Variance Std.Dev. Corr 
 Specific_habitat            (Intercept)        19.89788 4.4607        
                             log(Predator_mass)  0.23696 0.4868   -0.45
 Type_of_feeding_interaction (Intercept)        13.34598 3.6532        
                             log(Predator_mass)  0.08867 0.2978   -1.00
 Residual                                        3.95383 1.9884        
Number of obs: 34931, groups:  Specific_habitat, 16; Type_of_feeding_interaction, 5

Fixed effects:
                   Estimate Std. Error t value
(Intercept)         -3.6919     1.9849  -1.860
log(Predator_mass)   0.6254     0.1852   3.377

Correlation of Fixed Effects:
            (Intr)
lg(Prdtr_m) -0.755
convergence code: 0
Model failed to converge with max|grad| = 0.00435131 (tol = 0.002, component 1)

ranef(lm2)

$Specific_habitat
                                (Intercept) log(Predator_mass)
Coastal Bay                      -1.0293509         0.29755926
Coastal, SW & SE Greenland        6.1464920        -0.57594667
demersal food web                -0.2824898         0.29260144
estuary/coastal                   0.2565448         0.11056859
Euboikos and Pagassitikos Gulfs   5.4190217        -0.32491524
inshore                          -2.1204955        -0.19863649
nearshore                        -2.6864370        -0.10068252
Nearshore waters                 -0.1168369         0.24763310
Off the continental Shelf       -11.2551676         0.99021699
offShelf and on Shelf             5.0467912         0.73851030
open ocean                        6.1654872        -0.80127036
Pelagic                           0.3748829         0.06539027
Seasonal Pack Ice Zone           -4.4886450        -0.41345763
Shelf                             0.8120864        -0.15514253
Shelfbreak/open ocean            -0.4423849         0.11321355
transition region                -1.7994984        -0.28564208

$Type_of_feeding_interaction
                       (Intercept) log(Predator_mass)
insectivorous             3.400200         -0.2771532
piscivorous               1.356042         -0.1105321
planktivorous            -4.227637          0.3445982
predacious               -3.559123          0.2901071
predacious/piscivorous    3.030518         -0.2470201

tidy(lm2,conf.int = T)

fitted_lm2 <- augment(lm2) 

Deprecated: please use `purrr::possibly()` insteadDeprecated: please use `purrr::possibly()` insteadDeprecated: please use `purrr::possibly()` insteadDeprecated: please use `purrr::possibly()` insteadDeprecated: please use `purrr::possibly()` instead

ggplot(fitted_lm2, aes(x = log.Predator_mass., y = .fixed, colour = Specific_habitat)) +
    geom_line() 

Linear model 4

Assume different intercepts and slope for every feeding type.

lm4 <- lmer(log(Prey_mass) ~ 
                log(Predator_mass) + 
                (1 + log(Predator_mass)|Type_of_feeding_interaction), data = filtered_dataset)
summary(lm4)

Linear mixed model fit by REML ['lmerMod']
Formula: log(Prey_mass) ~ log(Predator_mass) + (1 + log(Predator_mass) |      Type_of_feeding_interaction)
   Data: filtered_dataset

REML criterion at convergence: 164001.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.5174 -0.7705  0.0110  0.6977  5.3503 

Random effects:
 Groups                      Name               Variance Std.Dev. Corr 
 Type_of_feeding_interaction (Intercept)        19.5231  4.4185        
                             log(Predator_mass)  0.1109  0.3331   -0.89
 Residual                                        6.3949  2.5288        
Number of obs: 34931, groups:  Type_of_feeding_interaction, 5

Fixed effects:
                   Estimate Std. Error t value
(Intercept)         -4.3564     1.9858  -2.194
log(Predator_mass)   0.6226     0.1548   4.022

Correlation of Fixed Effects:
            (Intr)
lg(Prdtr_m) -0.870

ranef(lm4)

$Type_of_feeding_interaction
                       (Intercept) log(Predator_mass)
insectivorous             1.487633        -0.15095072
piscivorous               4.435092        -0.42096904
planktivorous            -5.438602         0.25633769
predacious               -3.829539         0.37805913
predacious/piscivorous    3.345416        -0.06247707

tidy(lm4,conf.int = T)

Possible interaction effect for habitat and feeding type. In model 7, the slope is affected by the interaction btween habitat type and feeding type.

# #library(lme4)
# 
# #summary(lm7)
# ranef(lm7)
# tidy(lm7,conf.int = T)
# fitted_lm7 <- augment(lm7) 
# names(fitted_lm7)
#     ggplot(fitted_lm7, aes(x = log.Predator_mass., y = .fixed, colour = Specific_habitat)) +
#     geom_line() #+
#         #facet_grid(~ Specific_habitat)

AIC

library(MuMIn)

#model.sel(lm1, lm2, lm4, lm9, lm10, lm11,null, lm12, rank = AIC)

Take out model 5 – assumption doens’t make sense.

lm9 <- lmer(log(Prey_mass) ~ log(Predator_mass) + (1 | Type_of_feeding_interaction) + (1 | Type_of_feeding_interaction:Specific_habitat), data = filtered_dataset)
lm10 <- lmer(log(Prey_mass) ~ log(Predator_mass) + (1 + log(Predator_mass) | Type_of_feeding_interaction) + (1 + log(Predator_mass)| Type_of_feeding_interaction:Specific_habitat), data = filtered_dataset)

lm11 <- lmer(log(Prey_mass) ~ log(Predator_mass) + (1 + log(Predator_mass) | Specific_habitat) + (1 + log(Predator_mass)| Specific_habitat:Type_of_feeding_interaction), data = filtered_dataset)

# fitted_lm11 <- augment(lm11) 
# fitted_lm11
#     ggplot(nulllm11, aes(x = log.Predator_mass., y = .fitted)) +
#     geom_line(aes(colour = Specific_habitat))+facet_grid(~ Type_of_feeding_interaction) + null_line
# 
#MOVED THIS CODE TO THE BOTTOM BECAUSE NULLlm11 is at bottom

# fitdata11 <- ranef(lm11)[[1]] %>% 
#     rownames_to_column()%>% as_data_frame() %>% 
#     separate(rowname, into = c("habitat","feeding_type"), sep = ":") 
# colnames(fitdata11)[4] <- c("slope")
# fitdata11 %>% 
#     ggplot()+
#     geom_point(aes(x = feeding_type, y = fixef(lm11)[[2]] + slope,color = habitat) ) + 
#     geom_hline(yintercept = tidy(null)[2,2])

lm12 <- lmer(log(Prey_mass) ~ log(Predator_mass) + (1 + log(Predator_mass) | Specific_habitat) + (1 + log(Predator_mass)| Specific_habitat:Type_of_feeding_interaction) + (1 + log(Predator_mass)| Type_of_feeding_interaction), data = filtered_dataset)

# 
# fitdata12 <- ranef(lm12)[[1]] %>% 
#     rownames_to_column()%>% as_data_frame() %>% 
#     separate(rowname, into = c("habitat","feeding_type"), sep = ":") 
# colnames(fitdata12)[4] <- c("slope")
# fitdata12 %>% 
#     ggplot()+
#     geom_point(aes(x = feeding_type, y = fixef(lm12)[[2]] + slope,color = habitat) ) + 
#     geom_hline(yintercept = tidy(null)[2,2])

#model.sel(lm1, lm2, lm4, lm9, lm10, lm11,null, lm12, rank = AIC)

#anova(lm1, lm2, lm4, lm9, lm10, lm11,null, lm12)

#anova( lm9, lm10, lm11, lm12)

fitted_lm11 <- augment(lm11) 

Deprecated: please use `purrr::possibly()` insteadDeprecated: please use `purrr::possibly()` insteadDeprecated: please use `purrr::possibly()` insteadDeprecated: please use `purrr::possibly()` insteadDeprecated: please use `purrr::possibly()` instead

null <- glm(log(Prey_mass) ~ log(Predator_mass),data = filtered_dataset) 
fitted_null <- augment(null) %>% 
    rename(.fitted_null = .fitted)

Deprecated: please use `purrr::possibly()` insteadDeprecated: please use `purrr::possibly()` insteadDeprecated: please use `purrr::possibly()` insteadDeprecated: please use `purrr::possibly()` insteadDeprecated: please use `purrr::possibly()` instead

names(fitted_null)

[1] "log.Prey_mass."     "log.Predator_mass." ".fitted_null"       ".se.fit"            ".resid"             ".hat"              
[7] ".sigma"             ".cooksd"            ".std.resid"        

null_line <- geom_line(aes(x = log.Predator_mass., y = .fitted_null))
nulllm11 <- bind_cols(fitted_lm11, fitted_null)

fitted_lm11

    ggplot(nulllm11, aes(x = log.Predator_mass., y = .fitted)) +
    geom_line(aes(colour = Specific_habitat))+facet_grid(~ Type_of_feeding_interaction) + null_line

fitdata11 <- ranef(lm11)[[1]] %>% 
    rownames_to_column()%>% as_data_frame() %>% 
    separate(rowname, into = c("habitat","feeding_type"), sep = ":") 
colnames(fitdata11)[4] <- c("slope")
fitdata11 %>% 
    ggplot()+
    geom_point(aes(x = feeding_type, y = fixef(lm11)[[2]] + slope,color = habitat) ) + 
    geom_hline(yintercept = tidy(null)[2,2])

fitdata12 <- ranef(lm12)[[1]] %>% 
    rownames_to_column()%>% as_data_frame() %>% 
    separate(rowname, into = c("habitat","feeding_type"), sep = ":") 
colnames(fitdata12)[4] <- c("slope")
fitdata12 %>% 
    ggplot()+
    geom_point(aes(x = feeding_type, y = fixef(lm12)[[2]] + slope,color = habitat) ) + 
    geom_hline(yintercept = tidy(null)[2,2])

model.sel(lm1, lm2, lm4, lm9, lm10, lm11,null, lm12, rank = AIC)

Model selection table 
      (Int) log(Prd_mss)   class                                             random df    logLik      AIC    delta weight
lm12 -4.639       0.8521 lmerMod 1+l(P_m)|S_h+1+l(P_m)|S_h:T_o_f_i+1+l(P_m)|T_o_f_i 12 -72661.25 145346.5     0.00  0.740
lm11 -5.479       0.8311 lmerMod                  1+l(P_m)|S_h+1+l(P_m)|S_h:T_o_f_i  9 -72665.30 145348.6     2.10  0.259
lm10 -5.584       0.7093 lmerMod              1+l(P_m)|T_o_f_i+1+l(P_m)|T_o_f_i:S_h  9 -72670.70 145359.4    12.88  0.001
lm9  -5.532       0.7192 lmerMod                                T_o_f_i+T_o_f_i:S_h  5 -73533.87 147077.7  1731.23  0.000
lm2  -3.692       0.6254 lmerMod                      1+l(P_m)|S_h+1+l(P_m)|T_o_f_i  9 -73682.61 147383.2  2036.72  0.000
lm1  -5.689       0.9350 lmerMod                                       1+l(P_m)|S_h  6 -76536.60 153085.2  7738.69  0.000
lm4  -4.356       0.6226 lmerMod                                   1+l(P_m)|T_o_f_i  6 -82000.71 164013.4 18666.92  0.000
null -6.947       0.9966     glm                                                     3 -89077.81 178161.6 32815.11  0.000
Models ranked by AIC(x) 
Random terms: 
1+l(P_m)|S_h = ‘1 + log(Predator_mass) | Specific_habitat’
1+l(P_m)|S_h:T_o_f_i = ‘1 + log(Predator_mass) | Specific_habitat:Type_of_feeding_interaction’
1+l(P_m)|T_o_f_i = ‘1 + log(Predator_mass) | Type_of_feeding_interaction’
1+l(P_m)|T_o_f_i:S_h = ‘1 + log(Predator_mass) | Type_of_feeding_interaction:Specific_habitat’
T_o_f_i = ‘1 | Type_of_feeding_interaction’
T_o_f_i:S_h = ‘1 | Type_of_feeding_interaction:Specific_habitat’

anova(lm1, lm2, lm4, lm9, lm10, lm11,null, lm12)

refitting model(s) with ML (instead of REML)